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ABSTRACT 

Context, It is predicted that sources emitting UV radiation in the Lyman band during the epoch of reionization showed a series of 
discontinuities in their Lya flux radial profile as a consequence of the thickness of the Lyman line series in the primeval intergalactic 
medium. Through unsaturated Wouthuysen-Field coupling, these spherical discontinuities are also present in the 21 cm emission of 
the neutral IGM. 

^ — ' Aims. In this article, we study the effects these discontinuities have on the differential brightness temperature of the 21 cm signal 

| of neutral hydrogen in a realistic setting including all other sources of fluctuations. We focus on the early phases of the epoch of 

reionization, and we address the question of the detectability by the planned Square Kilometre Array. Such a detection would be of 
great interest, because these structures could provide an unambiguous diagnostic for the cosmological origin of the signal remaining 
after the foreground cleaning procedure. Also, they could be used as a new type of standard rulers. 

Methods. We determine the differential brightness temperature of the 21 cm signal in the presence of inhomogeneous Wouthuysen- 
Field effect using simulations which include (hydro)dynamics and both ionizing and Lyman lines 3D radiative transfer with the code 
LICORICE. We include radiative transfer for the higher-order Lyman-series lines and consider also the effect of backreaction from 
Q ' recoils and spin diffusivity on the Lya resonance. 

Results. We find that the Lyman horizons are difficult to indentify using the power spectrum of the 21 cm signal but are clearly visible 
on the maps and radial profiles around the first sources of our simulations, but for a limited time interval, typically Az as 2 at z ~ 13. 
Stacking the profiles of the different sources of the simulation at a given redshift results in extending this interval to Az as 4. When we 
take into account the implementation and design planned for the SKA (collecting area, sensitivity, resolution), we find that detection 
will be challenging. It may be possible with a 10 km diameter for the core, but will be difficult with the currently favored design of a 

5 km core. 

Key words. Radiative transfer - methods: numerical - intergalactic medium - large-scale structure of Universe - dark ages, reioniza- 
^ . tion, first stars 
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1 . Introduction Ho from Hubble Space Telescope observations, yield an optical 

depth t - 0.088 ± 0.014, which implie s a reionization redshift 
Cn Our universe was initially fully ionized, from the primordial very f z = i .6 ± 1.2 (iKomatsu et alJlfoHh . assuming an instanta- 
O ■ hot and dense phase until cosmological recombination, occur- neous transition . From these two observational results, we can 
ring at z ~ 1 100. After recombination, baryomc matter is made Mer that the EoR lasted mofjt nkely oyer an extended period . 
. of a neutral gas of hydrogen and helium . This period, known In addition to these two constraintS; new observations of the UV 
> ■ as the Dark Ages, lasted until the first light sources lit up at luminosity functions of galaxies in the reionization epoch have 
lower redshifts, as a consequence of the growth of primordial recently been published . These studies ^ in a position t0 put 
density fluctuations. These sources are responsible for the sub- constraints on the efficiency of galaxies in reionizing the inter- 
sequent reionization of the cosmological gas. Numerous ques- galactic medium (I GM). See for examp le the recent HUDF re- 
tions still have to be answered on that Epoch of Reionization suks at z ^ 7 _ 10 dBouwens et alJl2rMl2O10h . 
(EoR), and to date only a few observational results are in a po- 
sitio n to put constraints on t hat period. The Gunn-Peterson ef- The most promising probe of the EoR, however, is the fu- 
fect dGunn & Petersonll 19651) indicate s that the neutra l hydrogen ture observation of the 21 cm hyperfine line of neutral hydro- 
fraction is lower that 10" 4 at z < 5.5 dFan et aLl Hooefl Another gen (Hi). Indeed, Hi can be detected in emission or absorp- 
constraint is the optical depth due to Thomson scattering of cos- tion against the CMB at the wavelength of the redshifted 21 
mic microwave background (CMB) photons by free electrons, cm line corresponding to the transition be tween the two hyper- 
The seven-year results of the Wilkinson Microwave Anisotropy fine levels of the el ectronic ground state dHogan & Reeslll979| 
Probe (WMAP) experiment, combined with baryon acoustic os- IScott & Reeslll990l) . Pathfinder experiments such as LOFA 
cillation data from the Sloan Digital Sky Survey and priors on 
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1 For an interesting discussion about the redshift at which reioniza- 

tio n is complete, as deduced from the spectra of high redshift quasars, 

see Mesinaer (2010). 2 http://www.lofar.org/ 
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MWAS GMRlQ, 21CMA0 or PAPER0 will deduce information 
about the statistical properties of the signal, while the proposed 
Square Kilometre Array (SKAQ) will be able to deliver a full 
tomography of the IGM. 

An accurate prediction of the 21 cm signal requires the 
determination of the baryonic density, the ionization fraction, 
the kinetic temperature and the local Lya flux. Numerous pa- 
rameters can strongly affect the properties of t he predicted 
signal such as the size of the simulation box (ICiardi et al. | 
2003b] iBarkana & Loebl 120041 llliev et alJ 120061 iBaeketalj 
2009), the effect of the Lya wing scattering (ISemelin et al.l 



20071: IChuzhov & Zhend [2 007). or the nature of the reioniz 



ing sources (ICiardi et al.ll2003at [Mellema et alj|2006t llliev et alJ 



2007; IVolonteri & Gnedinl 120091: iBaek et al.ll2010h . Depending 
on these parameters, the EoR can extend over a variable red- 
shift range and produce a patchy morphology with sharper 
or smoother transition s between neutral and ionized periods 
dFurlanetto et al. 2006). Knowledge of the Lya flux is also im- 
portant, because it allows one to compute precisely the spin 
temperature of neutral hydrogen. Indeed, observations of the 
redshifted 21 cm line are only possible when the spin tem- 
perature is different from the CMB temperature, and among 
the physical processes likely to decouple these two tempera- 
tures, Jhe_ex2e£teddj}rmnan^ Wouthuysen-Field ef- 
fect dWouthuvsenl 119521: lFieldlll958l) . i.e. the pumping of the 
21 cm transition by Lya photons. For that reason, several 
studies focused on the Lya radiative transfer during the EoR 
(IBarkana & Loebl l200l purlanettol [2006L ISemelm et all l2007t 
IChuzhov & Zhendl2007tlBaek et al.ll2009h . 

In this paper we will consider the role of the radiative trans- 
fer of higher-order Lyman-series photons in the determination of 
the differential brightness temperature of the 21 cm signal. This 
is motivated by the fact that the local Lya photon intensity at a 
given redshift is made of two contributions: photons that were 
emitted below Ly/3 and redshifted to the local Lya wavelength; 
and photons emitted between Ly/3 and the Lyman limit, that red- 
shifted until they reached the nearest atomic level n. Due to the 
fact that the IGM is optically thick also to the higher levels, these 
photons were quickly absorbed and reemitted. But during these 
repeated scatterings the excited electrons pretty soon cascaded 
toward lower levels, a fraction of these radiative cascades end- 
ing with the emission of a Lya photon. 

The contribution of higher-order Lyman-line photons has 
been previously studied by several authors. However, while the 
first studies assumed a 100 perce nt efficiency for radiat i ve cas- 
cades to en d as a Lya photon ( e.g. Barkana & Loeb 2005), Hirata 
(120061) and lPritchard & Furlanettol(l2006l) showed that a sizeable 
fraction of the cascading photons ends in a two-photon decay 
from the 2s level. These authors calculated the exact cascade 
conversion probabilities and discussed the possibility of hyper- 
fine level mixing by scattering of Lyn photons. 

As a consequence of the thickness of the Lyman resonances 
in the primeval IGM, any primordial light source emitting pho- 
tons in the Lyman band is surrounded by discrete horizons, 
i.e. maximum distances photons with a given emission fre- 
quency can travel before redshifting into a Lyman resonance 
and being scattered by neutral hydrogen. Thus, Lya flux profiles 
around ionizing sources exhibit characteristic discontinuities at 
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the Lyman horizons. The importance of the Wouthuysen-Field 
effect in the determination of the differential brightness tempera- 
ture 6T\, of the 21 cm signal during the process of cosmic reion- 
ization motivates us to tackle the question of detecting corre- 
sponding brightness temperature discontinuities around the first 
light sources. In this paper, we include for the first time a correct 
treatment of the radiative transfer of higher-order Lyman-series 
lines into realistic numerical simulations of the EoR. The ques- 
tion is: will the horizons be observable in brightness temperature 
data over all the other sources of fluctuations, either physical 
(density, temperature, velocity effects) or instrumental? If this is 
the case, these spherical horizons of known radii at a given red- 
shift will provide a diagnostic to test systematic s and foreground 
removal procedures (see e.g. lHarker et al.ll2009l) . Without such a 
test, residual from imperfect removal may be impossible to dis- 
entangle from the cosmological signal. 

Our paper is organized as follows. In Sect. |2]we give some 
generalities about the Wouthuysen-Field effect and the differen- 
tial brightness temperature of the 21 cm transition. Lyman hori- 
zons are discussed in Sect. [3] Details on our numerical simula- 
tions are described in Sect. |4] while Sect. [5] shows our results. 
Finally, we set out our conclusions in Sect. [6] 

2. Wouthuysen-Field effect and differential 
brightness temperature of the 21 cm transition 

The differential brightness temperature of the 21 cm hyperfine 
transition is determined by the spin temperature T s , which is re- 
lated to the ratio between the densities of hydrogen atoms in 
the triplet («i) and singlet («o) hyperfine levels of the electronic 
ground state: 



"o go 



(1) 



where gi/go — 3 is the ratio of the statistical weights and 
T* = hvio/k^ - 0.0682 K is the temperature corresponding to 
the energy difference between the two hyperfine levels, with h 
the Planck constant, vio = 1420.4057 MHz the hyperfine transi- 
tion frequency and k% the Boltzmann constant. Three processes 
can excite these levels: absorption of CMB photons, collisions, 
and the Wouthuysen-Field effect, i.e. the mixing of the hyper- 
fine levels through absorption and spontaneous reemission of 
Lya photons. The latter effect dominates once the first luminous 
sourc es appear. Thus, the spin temperature can be written as (see 
e.g.[Furlanetto et al. 200(|: 



t- 1 + x a (Tfr l + x c r- 

1 I I 



(2) 



where T y is the CMB temperature, the effective color tem- 
perature of the UV radiation field, 7\ the gas kinetic temperature, 
x a the coupling coefficient for Lya pumping, and x c the coupling 
coefficient for collisions. 

The coefficient x c is given by: 



Al()Ty 



(C H + C p + C e ), 



(3) 



with Aio = 2.85 x 10 s the spontaneous emission factor of 
the 21 cm transition. Ce, C p and C e are the de-excitation rates 
due to collisions with neutral hydrogen atoms, protons and elec- 
tro ns respectively, for which we use the fitting formulae given 



trons respectivel y, tor w 
bv lKuhlenetai] (120061) . 
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The coupling coefficient for Lya pumping, x a , is given by: 



Xn — 



27 A w Ty 



(4) 



where P a is the Lya scattering rate per atom per second 
(Mad au et al.ll 1997b . In relation to the local Lye flux, the cou- 
pling coefficient can also be written as: 



167T 2 ?> 2 / t? 

21 A io TyirizC 



with e the electron charge, f a the Lya oscillator strength, »? e the 
electronic mass, c the speed of light and J a the angle-averaged 
specific intensity of Lya photons by photon number. The back- 
reaction factor S a accounts for spectral distortions near the Lya 
resonance, sourced by recoils and spin diffusivity. This factor 
cannot be neglected at the low kinetic temperatures characteriz- 
ing the very first stages of the EoR, when heating of the IGM by 
the first sources is still negligible. We include this important cor- 
rection factor using the simple analytical fit from Hirata (2006), 
accurate to within 1% in the range 7\ > 2K, T s > 2K, and 
10 5 < top < 10 7 , Tr,p being the Gunn-Peterson optical depth in 
Lya dGunn & Petersonll 19651) . 

When photons emitted with frequencies between Ly/3 and 
the Lyman limit redshift to the nearest atomic level n, they 
soon induce a radiative cascade, because of the thickness of 
the IGM to all Lyman levels. A fraction of the s e casc ades 
ends with the emission of a Lya photon. iHiratal (120061) and 
iPritchard & Furlanettol (120061) snowed that this process is not 
very efficient in producing Lya- photons, most of the radiative 
cascades ending in the 2s level, from which reaching the ground 
level is only possible through a two-photon decay. Furthermore, 
excitations of the 2s level via radiative excitations by CMB 
photons or collisional excitations are negligible at z < 400. 
As a consequence, these photons are lost for the Wouthuysen- 
Field effect. Besides, it should also be noted that direct contri- 
bution of Lyn scattering to the coupling of the spin tempera- 
ture to the gas temperature is not significant, as was shown by 
Pritchard & Furlanetto (2006), the number of scattering events 
before cascading being about 5, much less than the number of 
Lya scatterings through the line. 

Once the spin temperature is known, the differential bright- 
ness temperature of the 21 cm signal can be dete r mined 
by the following ex pression (see e.g. Ma dau et al.l Il997t 
Ciardi & Madau 2003): 



5T h = 28.1 x m (1 +6) 



Ob 



l+z\ 1/2 T s -T y 
10 



0.042 0.73 \ Q 



0.24 



1/2 



1 - 0.248 



mK, 



(6) 



where xui is the neutral hydrogen fraction, (1+6) the fractional 
baryon overdensity at redshift z, and Qb, £2 m , h, and Y p the usual 
cosmological parameters. Note that we do not consider in this 
paper any contribution from the proper velocity gradient. 

3. Lyman horizons 

Photons emitted by a source at redshift z s with frequency v s are 
observed with frequency v b s at z Q b s , with the following relation 
between emitted and observed frequencies: 



Vobs 



1 +Zs 
1 +Z b 



(7) 



1 000 
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Fig. 1. Evolution of the first five Lyman horizons, in comoving 
Mpc, as a function of the emission redshift. The horizon labeled 
Lyn is the maximal distance a photon emitted just below Ly(n+ 1) 
can travel. 



The comoving separation d between the source and the observer 
is determined by: 



d — c I df = c I 

J a® J H(z) 



dz, 



(8) 



with a the scale factor and H the Hubble parameter. At suffi- 
ciently high redshifts, the contribution of the cosmological con- 
stant to H can safely be neglected, and we can approximate 
H(z) ~ Hq(\ + z) VO m (l + z), with Ho the present-day value 
for the Hubble parameter. In that case, equation ((8} leads to: 



2c 



Hq yQ.„ 



:(1+Z.) 



-1/2 



(I \ -1/2 
/ Vobs ' 



1 



(9) 



The maximum distance a photon emitted just below Ly (n + 1 ) can 
travel before reaching the Lyn resonance is obtained when using 
v s = v n+ i and v b s = v n , with v„ the frequency corresponding to 
the transition from level n to the ground state: 



v„ = v LL (l - n 2 ), 



(10) 



where vll is the Lyman limit frequency. We then have the fol- 
lowing expression for the Lyn horizons: 



d„ = 



2c 



(l+Zs) 



-1/2 



Vn 

v„+i 



-1/2 \ 
- 1 



n > 2. 



(11) 



FigureQ]shows the evolution of the first five Lyman horizons, in 
comoving Mpc, as a function of the emission redshift. We see 
that photons emitted just below Ly/3 can travel more than 350 
comoving Mpc at z ~ 15, while photons emitted just below Ly^" 
cannot go farther than about 16 comoving Mpc. 

The existence of these Lyman horizons around primordial 
sources during the early EoR (while the Wouthuysen-Field cou- 
pling is not yet saturated) leads to discontinuities in their Lya 
flux profiles, which could create similar features in the bright- 
ness temperature profiles and maps. In order to answer quantita- 
tively this question, we have to consider a full modelling of the 
Lya pumping in a realistic numerical simulation of the EoR. 
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4. Numerical methods 

We simulate cosmological reionization as a three step process. 
First, a dynamical simulation is run using GADGET2 (Springel 
120051) . Snapshots of the simulation are produced at regular inter- 
vals. The interval is given in terms of the scale factor: Aa = 2~ 10 . 
Then the Monte-Ca r lo 3D radiative transfer code LIC ORICE 
(ISemelin et alJl2007b iBaek et alJl2009t Uliev et al.ll2009h is used 
to compute the ionizing continuum radiative transfer. We use 
two different dynamical simulations. The first one simulates a 
(100 h Mpc) 3 volume with 2 x 256 3 particles and assumes 
WMAP3 data alone cosmologi cal parameters: Q lrl = 0.24, Qa 
0.76. Q b = 0.042, h = 0.73 (ISpergel et al.ll2007l) . The source 
modelling (stellar sources, Salpeter IMF, with masses in the 
ra nge 1-120 M ) an d detailed reionization scheme, is described 
in lBaek et alJ (12010). where detailed analysis of the reionization 
history can be found. The second simulation is a box of size 
200 hT x Mpc, with 2 x 512 3 particles, such as the mass reso- 
lution is the same as in the first one. It assumes the more re- 
cent WMAP7+BAO+// cosmological pa rameters: fl m =0.272 , 
f2 A = 0.728, n h = 0.0455, h = 0.704 (iKomatsu et alJ[201lh . 
Finally, the Lyman-line radiative transfer part of LICORICE is 
used as a second post-process to study the Wouthuysen-Field 
effect, determine the spin temperature, and calculate the differ- 
ential brightness temperature of the 21 cm signal of neutral hy- 
drogen. In order to do so, we interpolated the output of the two 
reionization runs on a 256 3 (512 3 ) grid for the 100 h~ x (200 h~ l ) 
Mpc simulation respectively. In both cases, we emitted 1 .6 x 10 9 
Lyman band photons between each pair of snapshots. We em- 
phasize that the full radiative transfer in the lines is computed up 
to Lye. Using an isotropic r~ 2 profile around the sources would 
be much easier and faster. However, it has been shown that using 
the real, local Lya flux is important at early redshifts, for ioniza- 
tion fractions smaller than ~ 10% , when the Lya coupling is not 
yet full (IChuzhov & Zhenel2O07llBaek et al.ll2009l) . 

The implementation of radiative transfer for the higher- 
order Lyman-se ries fines is similar to the Lya radiative trans - 
fer described in ISemelin et all (120071) and IBaek et al.l (120091) . 
Important modifications are listed here: 

1. We now consider photons with frequencies in the range 
[v a , Vf], in order to include radiative transfer of the first five 
Lyman resonances (Lya - Lye). A 'flat' spectrum is assumed 
in this frequency range, i.e. every frequency has the same 
probability to be drawrfl We emit 1 .6 x 10 9 photons between 
each pair of snapshots. 

2. When calculating the optical depth of the propagating 
photon, the Lyn scattering cross-section is given by the 
Lorentzian profile: 



c«(v) = f u 

m e c 



(v-v») 2 + K/2) 2 ' 



(12) 



with fi„ the Lyn oscillator strength, the line center fre- 
quency, and Av£ the natural line width. Note that for photons 
whose frequency lies between Lyn and Ly(n + 1), we con- 
sider the possibility for scattering not only in the Lyn line 
wing, but also in the Ly(n +1) line wing. However, we find 



8 Given the narrowness of this frequency range, the use of a realis- 
tic non-flat spectrum results in minor differences, as we checked using 
blackbody spectra with thr ee different effective temp eratures: 5 x 10 4 
K, 10 s K, and 1.5 x 10 5 K. IChuzhov & Zhensl f2007) also showed that 
the spectral slope of the source does not have a strong influence on the 
slope of the scattering rate profile. 



that scatterings in the Ly(n + 1) line wing are rare and can 
safely be neglected. 

3. When a diffusion occurs in the Lyn (n > 2) line, we consider 
the probability for cascading to Lya. The probability for a 
Lyn photon to cascade is given by the value p = 1 - P„ p ^i s , 
where P np ^\ s is the simple diffusion probab ility from n to 
the ground state ( Pritchar d & Furla netto 2006). If the photon 
is not cascading, it is scattered by the atom. 

4. For cascading photons, we take into account the recycling 
fraction / leC ycie of radiative cascades that end in a Lya pho- 
ton emission. This fraction is zero for photons redshifting 
toward Ly/3 and about 1/3 for photons redshifting towar d 
Lyn, with n > 4 (lHiratall2006l: iPritchard & Furlanettoir2 006). 
The remaining radiative cascades do not contribute to the 
Wouthuysen-Field effect, since they end in the 2s level, from 
which they can reach the ground state only with a forbidden 
2y-decay. We compute scatterings until photons cascade or 
redshift within 2 thermal linewidths of a Lyn line center. At 
that distance from the center, given the thickness of the Lyn 
lines, photons are assumed to scatter without spatial diffu- 
sion and initiate a radiative cascade on the spot. 

5. Photons emitted below Ly/3 which reach the Lya resonance 
undergo top diffusions until they redshift out of the line, 
with: 



tgp 



3n H x m Al ya y 
2H ' 



(13) 



where ne is the proper number density of hydrogen, Aj_, ya the 
Lya wavelength, and y — 50 MHz the half width at half- 
maximum of the Lya resonance. Typical values are top ~ 
10 6 at z ~ 10. On the other hand, cascading photons are 
injected in the center of the Lya lin^| and thus undergo jTqp 
diffusions. 



In IBaek et al.l d2009l) . the photon propagation was stopped 
at 10 thermal linewidths of the Lya line center. Indeed, at that 
distance from the center, the photon mean free path is much less 
than 1 physical kpc, assuming standard density and temperature 
conditions for the redshifts of interest. However, the mean free 
path is about 4 physical kpc at 10 thermal linewidths of the Lye 
line center. We thus decided to stop the photon propagation at 2 
thermal linewidths, where the mean free paths are much shorter: 
~ 10~ 5 and ~ 7 x 10~ 4 physical kpc for Lya and Lye photons 
respectively. 

Note that we consider that Lya heating is neglig ible 
(IPritchard & Furlanettolf2 006: Furlanetto & Pritchardll2006b . In 



a recent paper, Meiksin (2010) studies heating by higher-order 
Lyman-series photons and gets heating rates as high as several 
tens of degrees per Gyr. However, we do not take this process 
into account, as our work concentrates on the very early phases 
of the reionization epoch. For the same reason, we completely 
neglect the effect of dust. While dust may play a significant role 
in Lyman-line radiative transfer during the later stages of the 
EoR, this is not expected in the beginning of the EoR and we 
consider a fully dust-free universe. 



9 Indeed, Furlanetto & Pritchard (2006) noted that even if the first ab- 
sorption is well blueward of the Lyman-series line center, the interme- 
diate states of the atom during the radiative cascade have small natural 
widths. 
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Fig. 2. Power spectrum of 6T\, at z = 11.05 for the fiducial 
200 h~ l Mpc simulation (solid line) and for a simulation neglect- 
ing higher-order Lyman-series radiative transfer (dotted line). 
The wavenumbers corresponding to the size of the three hori- 
zons are indicated. 



5. Results 

5.1. 3D power spectra 

The 3D power spectrum has been widely used as an efficient 
probe of the 5T\, fluctuations. However, we will show that this 
analysis tool cannot help us in detecting Lyman horizons in 
our simulated data. In Fig. [2] the 6T\, 3D power spectrum of 
our fiducial 200 h Mpc simulation (solid line) is compared, 
at z = 11.05, with the power spectrum of another simulation 
that does not take into account higher-order Lyman-series radia- 
tive transfer for the computation of 6Tb (dotted line). More pre- 
cisely, we propagate in the latter run the same number of pho- 
tons, but this time in the interval [v a , vp\. Arrows indicate the 
wavenumbers corresponding to the predicted sizes of the Lyy, 
LyS and Lye horizons, as given by equation (fTTT i. The redshift 
Z = 11.05 has been chosen because, as we will demonstrate 
later, this is when other analysis tools provide the best detection 
of the Lyman horizons. Also, the Hn bubbles at that moment in 
the simulations are numerous enough to have reliable statistics. 
As a consequence, the 3D power spectrum is expected to have a 
trustworthy shape and to be little affected by boxsize effects. On 
scales greater than the horizons, both spectra are very close to 
each other. The gain in power, compared with the fiducial sim- 
ulation, on scales shorter than the Lyman horizons, is caused by 
the combination of (i) the fact that, compared with the fiducial 
simulation, we do not loose any photons in the radiative cas- 
cades, and (ii) th e r~ 2 3 profile of P a on small scales due to Lya 
wing scatterings (Semelin et al. 2007; Chuz hov & Zh eng 2007). 
Note that the power spectrum show n in this work does not look 
like theoretical predictions (see e.g. lNaoz & Barkana 2008). We 
conclude from Fig.|2]that the two spectra have almost exactly the 
same shape, illustrating clearly the inability of 3D power spectra 
to help us in detecting Lyman horizons. 



5.2. Lyman horizons around individual sources 

Figure [3] shows the spherically averaged x a profile at z = 13.42, 
around the first source of the 100 h~ x Mpc simulation, appear- 




r [comoving Mpc] 



Fig. 3. Spherically averaged profile of the coupling coefficient 
for Lya pumping x a at z — 13.42, around the first light source 
in the 100 hr x Mpc simulation. The correct profile is repre- 
sented by the thick solid line, while the other two neglect back- 
reaction (dashed line) or radiative transfer for the higher-order 
Lyman-series lines (dotted line). As expected we observe steep 
decreases in the value of the profile at radii corresponding to the 
predicted Lyy, Ly<5 and Lye horizons (shown by arrows). 



ing at z — 14.06. The dotted line is the profile obtained when 
calculating radiative transfer including the Lya resonance only. 
In other words, this neglects the contribution of radiative cas- 
cades to the total Lya flux. The dashed line is the result of the 
correct radiative transfer calculation, including the five Lyman 
resonances from Lya to Lye, but without considering the effect 
of backreaction. Finally, the solid thick line is the correct aver- 
aged x a profile, with both Lyn radiative transfer and backreac- 
tion. The reason why there is no Ly/3 horizon is related to the se- 
lection rules forbidding Ly/3 photons to convert to Lya photons. 
We clearly observe discontinuities in the x a coefficient profile at 
the predicted positions. 

The corresponding profile for the differential brightness tem- 
perature is given in the upper panel of Fig. 2] using the same no- 
tation as in Fig. [3] Inside the tiny Hn bubble, i.e. in the first few 
comoving Mpc around the source (not shown in the figure), the 
signal is in emission. But outside this very small ionized region, 
the gas kinetic temperature is clearly lower than the CMB tem- 
perature at this early redshift. For that reason, the 21 cm signal 

appears in absorption. It reaches a minimal value of 190 mK 

at 4.3 comoving Mpc and is as large as about -60 mK at 10 co- 
moving Mpc from the source. Inclusion of higher-order Lyman- 
series radiative transfer leads to a significant increase by 25% of 
(—(57b) at 10 comoving Mpc compared to the case where only 
Lya radiative transfer is considered. Because the Wouthuysen- 
Field effect is dominant for determining the spin temperature, 
the discontinuities in the Lya flux translate into similar disconti- 
nuities in the differential brightness temperature profile. At this 
redshift, (—6T\,) decreases by about 5 mK, 2 mK and 0.2 mK at 
the Lye, Ly<5 and Lyy horizon locations respectively. This plot 
also shows the importance of taking into account backreaction 
for the correct computation of the spin temperature. Indeed, the 
box-average kinetic temperature is (TV) = 4.63 K at z - 13.42, 
and at such a low temperature the backreaction correction factor 
is about 0.6. In the lower panel of Fig. [4] we plot the gradient 
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Fig. 4. Upper panel: Spherically averaged profile of the differ- 
ential brightness temperature at z — 13.42, around the first light 
source in the 100 Mpc simulation. The same notation is used 
as in Fig. [3] Note that we plot the opposite value of 6T\, . Lower 
panel: Gradient of the spherically averaged profile. 



of the spherically averaged profile. Steep decreases in the tem- 
perature profile at the positions of the Lyman horizons result in 
prominent peaks in its gradient that make these horizons even 
easier to detect. The most important conclusion from Fig. |4]is 
that the different steps in the x a profile at the Lyman horizons 
correspond to similar, clearly seen decreases in the differential 
brightness temperature profile. 

We show in Fig. [5] a map of the quantity -5T b x r 2 around the 
same source at the same redshift z = 13.42, with r the distance to 
the source center. Mapping this quantity enhances visualization. 
Indeed, from equations d2]i and ©, and since[3 x a oc r~ 2 , it can 
easily be shown that 5T\, is also proportional to r~ 2 . Mapping the 
product -6T\, x r 2 will thus straighten up the radial profiles, mak- 
ing the horizons more apparent. The slice thickness is 2 bT l Mpc 
and the map scale is logarithmic. The three spherical, concentric 
horizons are marked with white, yellow and red arrows for the 
Lye, Lyc) and Lyy discontinuities respectively. On this map, the 




10 More precisely, taking wing scatterings into a ccount slightly steep- 
ens the profile on small scales, as noted above ( Semeli n et al.l [20071 : 
IChuzhov & Zhendl2007l) . 



Fig. 5. Map of the quantity -6T b x r 2 at z = 13.42, with r the 
distance to the source center, in arbitrary units. The Lye, Ly<5 
and Lyy horizons are marked with white, yellow and red arrows 
respectively. The size of the box is 100 h _i Mpc and the slice 
thickness is 2 hr l Mpc. The color scale is logarithmic. 



Lyt5 horizon is particularly obvious. Note the perfectly spherical 
shape that characterizes these features. 

5.3. Effect of the nearby sources 

From Figs. [3] to [5] we are able to conclude that the existence 
of Lyman horizons around the first light sources create similar 
features in the 21 cm signal at the very beginning of the EoR. 
However, these weak discontinuities, whose amplitude is only 
a few mK, will be progressively affected by the contribution of 
other nearby sources located closer than the ~ 50 comoving Mpc 
corresponding to the farthest Lyy horizon. Also, Lyman photons 
emitted between Lya and Lyy by more distant sources can travel 
more than 50 comoving Mpc, and thus are likely to reach the 
outskirts of many other sources. The increasing influence of the 
other UV emitters on the brightness temperature profile of an 
individual source is shown in Fig.[6]for the 100 h Mpc simu- 
lation box, from z = 13.22 to z = 11.64. This redshift interval 
is covered by 10 snapshots and the number of emitting cells is 
6 at z = 13.22, 15 at z = 12.66, 38 at z = 12.13, and 86 at 
Z = 1 1 .64. Early in the simulation, the three Lyman horizons are 
clearly seen as prominent peaks in the gradient of the tempera- 
ture profile (upper left panel). Later, as the influence of the other 
sources becomes stronger, fluctuations in the profile grow until 
they reach the same amplitude as the horizons. At the beginning 
of this process, some of the horizon peaks are still observable 
(upper right and lower left panels), but after a certain time they 
are not distinguishable anymore (lower right panel). We find that 
for this particular simulation the Lye and Ly<5 horizons are ob- 
servable during a redshift interval Az ~ 2 and that the size of 
these discontinuities is about 2-4 mK. The Lyy horizon is visi- 
ble for a shorter period Az < 1.5. 
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Fig. 6. Gradient of the spherically averaged profiles of the differential brightness temperature, around the first light source appearing 
at z — 14.06 in the 100 hr l Mpc simulation. The four panels show the profile at different redshifts, from z = 13.22 to z — 11 .64. 
At each redshift, arrows show the predicted position of the Lye, Ly6 and Lyy horizons. The Lye and Ly6 horizons can be detected 
during a redshift interval Az ~ 2. The fainter Lyy horizon is visible for a shorter period Az < 1 .5. 



The size of our simulation box can also be important in 
this attenuation process. Because of the periodic boundary 
conditions, a photon emitted, say, between Lya and Ly/5 can 
travel a distance greater than the box size, resulting in a self- 
contamination of the Lya flux. Also, the smaller the simulation 
box, the more homogeneous is the distribution of sources around 
any given point in the box. In order to reduce these effects, we 
considered the larger, 200 h Mpc simulation and find that, 
despite the larger volume, the Lyman horizons around the first 
source of the simulation are observable during about the same 
redshift range as for the first source in the 100 h~ l Mpc simula- 
tion: Az ~ 2 for the Lye and Lyd discontinuities, Az ~ 1.5 for 
the Lyy step. These are maximal intervals and sources lighting 
up later will have their horizons observable for a reduced period. 
As an example, the seventh source in the 200 h Mpc simula- 
tion, appearing at z = 13.03, has Az ~ 1.2—1.5 for Lye and Ly<S, 
and the Lyy horizon can be detected in the first few snapshots 
only. 

The Lye and Ly<5 horizons are thus the best candidates for 
detection, the latter one being visible a little bit longer. The Lyy 
discontinuity is weaker because of the r~ 2 scaling of the Lya 
flux, and is thus more sensitive to the influence of the other 
sources. 



The growing fluctuations appearing in the spherically aver- 
aged temperature profiles clearly originate from the Lya back- 
ground that develops as more and more sources light up. At first 
sight, density or gas temperature fluctuations could also be sus- 
pected (see equation|6]l. We thus ran test simulations in which we 
computed 6T\, assuming a perfectly homogeneous and isother- 
mal IGM in the Lyman radiative transfer post-processing phase. 
The density was chosen to be the critical density and the gas ki- 
netic temperature was assumed to be equal to the box-averaged 
value of the normal run. In this way, only the Lya flux was likely 
to let 6T\, vary from one point to the next. The results of this test 
showed that the Lyman discontinuities are not seen for a signif- 
icantly longer period than in the standard case. Indeed, we were 
able to identify the three horizons in a few more snapshots only, 
and we conclude that the density and temperature anisotropics 
around the source are not responsible for the disappearance of 
the Lyman features. A final possibility in explaining the vanish- 
ing of the Lyman steps is Monte-Carlo noise. Indeed, the total 
number of photons per snapshot is constant, but the source num- 
ber grows quickly. As a result, the number of photons per source 
decreases very rapidly. However, using very expensive simula- 
tions with ten times more photons between each pair of snap- 
shots (16 x 10 9 ) in both the 100 ft -1 and 200 ft -1 Mpc simula- 
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Fig. 7. Gradient of the stacked ST-o radial profile at z — 1 1 .05 
for the 200 h Mpc simulation. The stacked profile is obtained 
by averaging the individual profiles of all the source-containing 
cells at this redshift. This method allows us to detect the Lye and 
Ly<S horizons, whose predicted positions are indicated by arrows. 



tions results in undistinguishable radial profiles. We thus clearly 
demonstrated that the Lyman horizons disappear around individ- 
ual sources because of the growing influence of the Ly<? flux of 
the nearby sources. 

5.4. Stacking 

In order to extend the detectability of the Lyman steps to lower 
redshifts, we considered a stacking technique. This allows one 
to smooth the fluctuations present in the individual profiles and 
to strengthen the visibility of the horizon steps. As an example 
to illustrate the efficiency of this method, we stacked the 5T\, 
radial profiles around all the cells that contain sources in the in- 
terpolation grid of the 200 h Mpc simulation, at z — 11.05, 
when the Lyman horizons are undetectable around the vast ma- 
jority of the individual sources. For some sources, small features 
are seen at the predicted horizon positions, but are indiscernible 
from fluctuations of similar amplitude. At this redshift we find 
1433 source-containing cells. The individual profiles are sim- 
ply added and the result is divided by the number of cells we are 
considering! FigureQshows the gradient of the stacked profile. 
The presence of the Lye and Lyd> horizons is now obvious. We 
find that this method allows us to detect these two discontinuities 
unambiguously down to redshifts close to 10. Interestingly how- 
ever, the stacking method does not provide good results for the 
smaller simulation. We believe that this comes from the number 
of source-containing cells being about five times smaller in the 
100 hT x Mpc run at a given redshift. This has the consequence 
that the number of profiles that are averaged is not enough to 
smooth the individual fluctuations. 



11 We can safely ignore the fact that the stacked sources have different 
peculiar velocities perpendicular to the line of sight. Indeed, a rapid 
estimation shows that a peculiar velocity about 300 km s -1 leads to a 
redshift difference Az = (1 + z) v/c ~ 10~ 2 . As shown in Fig. [TJand 
equation dl lb . this results in completely negligible differences in the 
horizon sizes. 



Fig. 8. Gradient of the 5T\, profile around the first source appear- 
ing in the 200 h~ x Mpc simulation at z — 12.84, with (dotted 
line) and without (solid line) SKA-like noise. 

5.5. Detectability 

We would now like to examine the possibility for detecting the 
Lyman horizons with the planned Square Kilometre Array. To 
reach this goal we have to include in our simulated data both 
the instrument resolution and noise. We choose the nois e power 
spectrum given by equation (12) of ISantos et al.l (1201 ll) . which 
models the resulting noise after foreground removal has been 
performed: 

, nA 2 D 2 m . d J 2 

P N (k,0) = r 2 y %LJ*L, (14) 

f 0^tot 

where k is the wavenumber, 9 the angle between the wave vec- 
tor and the line of sight, r{z) the comoving distance to redshift 
z,y a conversion factor between frequency intervals and comov- 
ing distances, A the observation wavelength, Z) max the maximum 
baseline, r sys the system temperature, fo the total observation 
time, and A tot the total collecting area. Following Santos et al., 
we use D max = 10 km, fo = 1000 hours and a system tempera- 
ture modelled by their equation (13) as the sum of the receiver 
noise temperature (assumed to be 50 K) and the sky temperature 
(dominated by the Galactic synchrotron emission): 

/ 1 \ 2 ' 55 

r sys (z)=50 + 60l^!j K. (15) 

We compute the total collecting area using the sensitivity value 
of 4000 m 2 /K, in agreement with the design reference for the 
instrument. Following these authors, we also use 70% of A tot in 
the noise calculation, considering that the rest of the collecting 
area is used for point source removal and calibration. 

Figure|8]illustrates the dramatic effect that noise may have on 
the detectability of Lyman horizons. The plotted example shows 
the gradient of the 6T\, radial profile around the first source of 
the 200 h~ x Mpc simulation at z = 12.84, with (dotted line) and 
without (solid line) noise. While the Lye and Ly£ horizons are 
clearly detected in the noise free profile, this is no longer true 
when noise is added. 

However, the field of view (FoV) subtended by our simula- 
tion (less than 3 deg 2 in the range z — 1 1 - 14) is much smaller 
than probable SKA field of view. If the ratio between both fields 
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is N, then noise will be reduced in the stacking procedure by a 
factor y/N larger for the SKA FoV than for the simulation FoV. 
We will here consider a field of view of 400 deg 2 for the SKA. 
To simulate the efficiency of the stacking procedure for such a 
large field of view, we decrease the noise level in the simula- 
tion FoV by the appropriate y/N factor. The fact that we stack 
over the volume of our simulation cube is not problematic in this 
analysis. Indeed, a short calculation shows that the full length 
of our 200 h' 1 Mpc simulation represents a redshift difference 
Az ~ 1.5, which corresponds to a ~ 6% difference in the horizon 
sizes at z ~ 13, according to equation ( fTTT i. For example, this 
represents 1 .6 comoving Mpc for LyS, or about 0.5'. 

The most limiting factor will probably be the resolution of 
the instrument. The proposed SKADS-SKA implementation for 
observations between 70 MHz and 450 MHz suggests a 5km- 
diameter cor^El This gives a resolution larger than 2' in the red- 
shift range z — 1 1 - 14 relevant for our study. In order to mimic 
the effect of the limited instrumental resolution, we apply 3EB 
gaussian smoothing with a FWHM equal to the resolution. We 
compare in Fig. [9] a map at z = 11.05 of the full-resolution, 
noise-free simulated data (left panel) with the same field when 
one includes noise and convolution with the appropriate resolu- 
tion (right panel). The full-resolution map shows a characteristic 
signature of the early EoR 21 cm signal in the vicinity of the 
sources, with practically spherical Hn bubbles being surrounded 
by absorption regions with deep 6Tb local minima. On the other 
hand, only some of the largest bubbles are still detected in the 
convolved map. This will have an influence on the identification 
and localization of the first sources, and hence on attempts to ob- 
serve Lyman horizons. Also, the limited resolution of the future 
observations will cause a blurring of the small discontinuities we 
want to detect. 

Finding an algorithm to efficiently detect the Hn bubbles in 
the convolved maps will be the subject of a future paper. Here, 
we tackle the question of pinpointing the source positions on the 
sky beyond the resolution accuracy, which could leave its mark 
on the stacking procedure, in the following way. Instead of using 
the exact source position when building the stacked profile, we 
randomly distribute its location inside a sphere centered on the 
exact position. The diametre of that sphere corresponds to the 
SKA resolution at the considered redshift. Figure [TOl shows the 
gradient of the stacked profile assuming the observation strategy 
described above (thick solid line), compared to the case without 
noise and with full resolution, and in which the position of the 
sources is perfectly known (dotted line). The important field of 
view allows us to significantly reduce the effect of the noise, but 
the 2' resolution has a dramatic effect on the Lyman horizons, 
which are completely wiped off. The dashed line shows that with 
a 1' resolution (10 km core), the Lyd horizon is detected as a 
wide and weak hump. In order to maximize the detection prob- 
ability, we assumed here perfect foreground removal and chose 
a redshift (z = 13.42) at which the Lyman horizons are not yet 
contaminated by the other nearby sources. We deduce from this 
figure that detection of the Lyman horizons in the early phases of 
the EoR using SKA observations will probably be possible only 
with a 1' resolution. 



12 http : //www . skads-eu . org/PDF/SKADS Jfhite_Paper_lQ83 18 
_dist.pdf 

13 For simplicity reason, we consider the same resolution in the fre- 
quency direction as in the other two spatial directions. 




r [comoving Mpc] 

Fig. 10. The thick solid (thin dashed) line shows the gradient 
of the 6T\, stacked profile in the 200 hr x Mpc simulation at 
z — 13.42, with SKA-like noise, assuming a field of view of 
400 deg 2 , and a 2'30" (1' 15") resolution. For comparison, we 
also show the gradient of the 6T\, stacked profile, at the same 
redshift, without noise nor resolution effect (dotted line). 

6. Discussion and conclusions 

In this paper we studied the effects of higher-order Lyman-series 
radiative transfer on the differential brightness temperature of 
the 21 cm signal of neutral hydrogen during the early stage of 
the EoR. This signal was computed in the presence of inhomo- 
geneous Wouthuysen-Field effect using the Monte-Carlo 3D ra- 
diative transfer code LICORICE. We showed that the disconti- 
nuities in the Lye flux radial profile of the first light sources, 
caused by the existence of discrete Lyman horizons, lead to sim- 
ilar features in the differential brightness temperature profile. 

We found that the Lye and Ly<5 horizons are detected in our 
simulations around individual sources during a redshift interval 
Az ~ 2 after the first source lights up. Later, the Lya flux at a 
given point has a growing smooth background component due to 
the cumulative contribution of the other sources, and the Lyman 
steps fade away. However, we showed that stacking the indi- 
vidual profiles makes the visibility period significantly longer 
(As- 4). 

An important caveat has to be mentioned here. We have con- 
sidered in our simulations purely stellar sources. However, it 
should be kept in mind that considering an important contribu- 
tion from quasars would result in the Wouthuysen-Field effect 
being significantly induced by Lya- photons coming from colli- 
sional excitation of hydrogen atoms by secondary electrons pro- 
duced by X-rays. If an important part of the total number of Lya 
photons is not coming from redshifting and cascading photons, 
the interesting features in 6T\, we discussed in this work would 
be more difficult to extract. 

On the other hand, let us note that the size of the Lyman dis- 
continuities could increase if the number of Ly a photons from 
radiative cascades is enhanced. In a recent paper, Meiksin (2010) 
claims that the number of such photons could be up to 30% 
higher than previously considered. In this case, the size of the 
steps in the x a profile would increase interestingly, together with 
the size of the corresponding 6T\, discontinuities. 

Detections of these discontinuities would be of first impor- 
tance for two reasons. First, such features would be very inter- 



9 



P. Vonlanthen et al.: Distinctive rings in the 21 cm signal of the epoch of reionization 




Fig. 9. Maps of 6T\, at z = 11.05. The slice is 200 h~ l Mpc on a side and has a thickness of 2 Mpc. The scale, in units of mK, 
is linear. The left panel shows our simulated data cube at full resolution. The right panel shows the effect of instrumental noise and 
instrumental resolution. At that redshift, the latter is slightly larger than 2'. 



esting for researchers involved in foreground removal. Indeed, 
removing the much brighter foregrounds (Galactic synchrotron 
and free-free emission, extragalactic sources, ionospheric distor- 
tions) will be a serious issue for future 21 cm observations. But 
Lyman horizons are structures whose size is precisely known 
at a given redshift (or frequency), and for that reason posi- 
tive detections would ensure that the laborious removal task 
has been correctly handled. The second reason is the possi- 
ble use of these horizons as standard rulers, as was noted by 
iPritchard & Furlanettol(l2006l) . 

When assessing the question of detecting Lyman horizons 
with the planned Square Kilometre Array by taking into account 
the reference implementation and design, we found that such a 
detection will be very challenging. Indeed, the horizons are not 
detected with a 2' resolution. A 1' resolution may allow us to 
detect the Ly£ horizon. 
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